Contractor-Renormalization approach to frustrated magnets in a magnetic field 
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We propose to use the Contractor Renormalization (CORE) technique in order to derive effective 
models for quantum magnets in a magnetic field. CORE is a powerful non-perturbative technique 
that can reduce the complexity of a given microscopic model by focusing on the low-energy part. 
We provide a detailed analysis of frustrated spin ladders which have been widely studied in the 
past: in particular, we discuss how to choose the building block and emphasize the use of their 
reduced density matrix. With a good choice of basis, CORE is able to reproduce the existence or 
not of magnetization plateaux in the whole phase diagram contrary to usual perturbation theory. 
We also address the issue of plateau formation in two-dimensional bilayers and point out the analogy 
between non- frustrated strongly anisotropic models and frustrated SU(2) ones. 

PACS numbers: 75.10.Jm, 75.60.Ej 
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I. INTRODUCTION 

In the presence of a magnetic field, quantum magnets 
exhibit fascinating properties. In particular, it can hap- 
pen that the uniform magnetization along the field ex- 
hibits plateaux for rational values, which has given rise to 
lots of theoretica l 1 ' 2 ' 3 ' 4 and experimental^^ work. The 
appearance of such plateaux was found to be favored by 
magnetic frustration. 

More recently, experiments on spin dimer compounds 
(such as TlCuCl 3 , KCuCl 3 and BaCuSi 2 6 , see Ref. i, 
Ifl flolfrTh have shown that the triplet excitations behave 
as bosons that can form superfluid (absence of plateau) 
or crystalline (finite plateau) phases depending on the 
competition between repulsive and kinetic interactions. 
Moreover, the possibility of having both orders, namely a 
supersolid, could potentially be observed in related com- 
pounds. As is well-known, frustration reduces triplet de- 
localization and thus, is favorable to solid behaviour, i.e. 
plateau formation, or supersolid behaviour. However, 
frustrated spin models are difficult to study numerically 
due to the sign-problem of the Quantum Monte-Carlo 
(QMC) technique and the absence of reliable large-scale 
numerical techniques in two dimensions or higher. How- 
ever, the effective bosonic models themselves can often 
be simulated when the frustration has disappeared and 
is absorbed in the effective parameters; indeed, such ef- 
fective bosonic models have been proposed either based 
on perturbation theor y 2 i 12 i 13 or on phenomenological 
grounds^ Therefore, we think that it would be highly 
desirable to derive non-perturbative effective parameters 
directly from microscopic models and we propose to use 
the contractor-renormalization (CORE) technique to do 
so. 

The CORE method has been proposed in Ref.[lll6|Il3 
as a systematic algorithm to derive effective Hamiltoni- 
ans and operators that contain all the low-energy physics. 
In principle, these effective operators are given by an 
infinite cluster expansion. CORE has been successfully 
applied to a variety of both magneti o 17 i 18 ' 19 i 20 ' 21 i 22 and 
dope d 23 i 24 low-dimensional systems and it turns out that, 



in most cases, the cluster expansion converges quite fast, 
which is a necessary condition for any practical imple- 
mentation of this algorithm. Still, the issue about CORE 
convergence is crucial and currently debated* 2 ^ 

In section HH we remind the reader of the CORE algo- 
rithm and investigate the frustrated 2-leg ladder. Being a 
well-known model, we can compare our findings to other 
well-established numerical results and discuss the accu- 
racy of CORE as well as its convergence. In section Hill 
we investigate how to choose the best block decomposi- 
tion and how to select the low-lying states to keep by us- 
ing information obtained with the exact reduced density- 
matrix of a block embedded in a large system. Finally, in 
section ITVl we turn to some two-dimensional (2D) bilayer 
spin models which are candidates for observing some of 
these exotic bosonic phases. 



II. FRUSTRATED 2-LEG LADDER 



The Hamiltonian Ti of the spin S = i frustrated anti- 
ferromagnetic Heisenberg ladder in an external magnetic 
field h reads: 



H = J± ^ S r 
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In accordance with Fig. [TJ the index r — 1, . . . , L repre- 
sents the L different rungs of the ladder whereas i = 1, 2 
indicates the two chains which constitute the ladder and 
we use periodic boundary conditions along the legs. The 
rung spin exchange J± is set to 1, while J x and Jd stand 
for the interactions along the chain and diagonal respec- 
tively, which makes the ladder a frustrated system. Note 
that due to symmetry J x and Jd can be interchanged. 
The role of frustration in the plateau formation can also 
be understood in a related ladder models 

Throughout this paper, we only consider the physical 
properties at zero temperature and in the presence of a 
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finite magnetic field, or more specifically, the possibility 
or not of a magnetization plateau at half-saturation : 
m z = \m sat . 



A. Summary 

The frustrated Heisenberg ladder has been studied 
with various analytical and numerical technique o 27 ' 28 ' 29 . 
In the absence of magnetic field, there are two main re- 
gions called rung singlet phase and Haldane phase. The 
effect of a magnetic field has been studied with an ex- 
act diagonalization (ED) technique 30 in order to clarify 
the presence or not of finite-magnetization plateaux. For 
all the Hamiltonians that we will consider (including the 
effective ones), S 1 ' ' is a good quantum number; there- 
fore, it is sufficient to compute the ground-state energy 
in all S z sectors (in the absence of any magnetic field) 
and then perform a Legendre transform to get the full 
magnetization vs field curve m z {h). We have applied the 
same ED technique and we provide on Fig. [5] a sketch 
of its phase diagram. Data have been extrapolated to 
the thermodynamic limit after standard finite-size scal- 
ing analysis of the plateau size (see examples on Fig. [3]). 
A large magnetization plateau phase is found around the 
strongly frustrated region J x ~ Jd- We note the exis- 
tence of different phases without plateau; in particular 
at large Jd ~ J x , there is a first order transition to the 
so-called Haldane phase. 

For the plateau phase, we draw on Fig.[2]a typical mag- 
netization curve obtained on a 2 x 16 ladder. It exhibits 
singularities at critical fields 3 -^ and a large half-saturation 
plateau. Note that in order to mimic the thermodynamic 
limit, we have drawn a line connecting the middles of the 
finite-size plateaus of the finite system. 



B. Perturbation theory 



FIG. 2: (Color online) Size of the magnetization plateau at 
m z = \rn aa t- Results have been obtained after finite size 
scaling analysis of exact diagonalization data obtained on 2 x 
L ladder (up to L = 20). 
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FIG. 3: (Color online) Finite-size scaling of the size of the 
magnetization plateau at m z — ^m sa t- Data from numer- 
ical exact diagonalization is shown for an example of a fi- 
nite plateau (J x — 0.7, Jd = 0.5) and a vanishing plateau 
(Jx = 0.4, J d = 0.1). 



a perturbation theory^ 2 -. Using pseudo-spin S = \ oper- 
ators a r for these two states, Eq. (fTJ) can be rewritten 
as an effective Hamiltonian H e ff. Proceeding further, a 
Jordan- Wigner transformation leads to a system of inter- 
acting, one-dimensional (ID), spinlcss fcrmions: 



When the only nonzero coupling is Jj_, the ground- 
state of the ladder is simply the product of singlets on 
each rung. The states of a given rung are labelled as a sin- 
glet \s) r = 7^(1 Ti)r — I IT)r) an d the three components 
of a triplet : |^) r = | U) r , \t ) r = ^(| U)r + I IT)r) 
and = | TT)r- I n the presence of a finite mag- 

netic field, due to Zeeman splitting, one can restrict the 
Hilbert space on each rung to |s) and and then do 



H t v = t y~] (c|,c r+ i + h.c.) + V n r n r+1 



L 



(2) 



where t describes the hopping, V the nearest neighbor 
interaction and [i the chemical potential which can all be 
expressed in terms of the previously introduced interac- 
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FIG. 4: (Color online) Magnetization m z along the field (nor- 
malized to its saturation value m S at) as a function of magnetic 
field h on a 2 x 16 ladder with J x = 0.7, and Jd — 0.5; data 
from numerical exact diagonalization has been used. 

tions J_l, J x , and Jd : 

, Jx Jd t r Jx "T" Jd /„x 

« - - w 

[X = J± — h 

In the particle language, the occurence or not of a plateau 
translates into the existence of single-particle gap and the 
magnetic field plays the role of an effective chemical po- 
tential /i (see Ref.y). Since the metal-insulator transition 
of the t — V model 3 -^ occurs at half- filling when V/\t\ = 2, 
we conclude that there will be a finite plateau at half its 
saturation value when 

i < Jd/Jx < 3. (4) 

This property is shown in Fig. [5] and reasonably agrees 
with the exact results of Fig. [21 although we note quanti- 
tative differences in the non-perturbative regime. There- 
fore, we now turn to a non-pertubative CORE approach 
with the same basis. 



C. CORE approach in the rung basis 

The Contractor Renormalization (CORE) method has 
been formulated in 1994 by Morningstar and Wein- 
stei n 15 i 16 and has been used subsequently to study both 
magnetic system o 17 i 18 i 19 i 20 i 21 i 22 and doped one o 23 i 24 ' 33 . 
The idea of this non-perturbative method is to derive an 
effective Hamiltonian within a truncated basis set which 
allows to reproduce the low energy spectrum. This means 
that the original model is replaced by a model with fewer 
states but a more complicated Hamiltonian under the 
condition that the retained states of the modified model 
have an overlap with the set of lowest lying eigenstates 
of the full original theory. 



FIG. 5: (Color online) Phase diagram using either perturba- 
tion theory or range-2 CORE. Reduced density matrix calcu- 
lations on a 2 x 12 ladder indicate that, for large J x and Jd, 
the reduced density matrix weight of the singlet becomes very 
small while the S z = triplet weight increases substantially 
(region above the blue dashdotted line, see text for details). 
For comparison with the exact results, see Fig. [2] 

For clarity, we briefly remind the main CORE steps 
and refer to the literature for more details. First, one 
needs to choose a basic cluster and diagonalize it. Then 
M low-energy states are kept and the remaining states 
are discarded. Generally, the M lowest states are re- 
tained, but this is not a necessity as we will discuss in 
Sec. IIIII The second CORE step is to diagonalize the full 
Hamiltonian H on a connected graph consisting of N c 
clusters and obtain its low-energy states \n) with ener- 
gies e n . Thirdly, the eigenstates |n) are projected on the 
tensor product space of the retained states and Gram- 
Schmidt orthonormalized in order to get a basis |\l/ n ) of 
dimension M Nc . Fourthly, the effective Hamiltonian for 
this graph is defined as 

h Nc = e »l*n>(*n|. (5) 

n=l 

Fifthly, the connected range-iV c interactions h^ nn can 
be calculated by substracting the contributions of all con- 
nected subclusters. And finally, the effective Hamiltonian 
is given as a cluster expansion as 

h core = j2 hl + j2 +Y, h ^ n + ---- ( 6 ) 

i <ij> <ijk> 

Then, of course, one has to study this new effective 
model, which can still be a difficult task. One possi- 
bility is to iterate CORE in the renormalization group 
spirit and study the properties of its fixed point. Here, 
following Ref. [20| . we propose to check the validity of 
the effective Hamiltonian after one step by comparing its 
properties to the exact ones on a given finite cluster. Of 
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FIG. 6: (Color online) For J x = 0.4, J d = 0.2 the norm of 
hf nn quickly decreases with increasing range r, whereas for 
J x = 0.8, Jd = 0.7 the norm of h c j onn has the same order of 
magnitude for ranges 2 < r < 8. 



FIG. 7: (Color online) With increasing range, the CORE re- 
sults converge towards the exact result. Parameters are: 2xL 
with L — 14; J x = 0.4 and Jd = 0.2. The effective model was 
obtained up to range 5, and then solved on a 14-site chain. 



course, once the effective Hamiltonian has been shown to 
be accurate, it can be simulated exactly on much larger 
lattices than the original model, or thanks to other nu- 
merical techniques. For instance, although the original 
model is frustrated and cannot be simulated efficiently 
by Quantum Monte-Carlo (QMC), there are cases where 
the effective Hamiltonian can. 

We begin our CORE considerations for the spin S = 
i antiferromagnetic Heisenberg ladder with a rung as 
simplest possible basic cluster. Because of the Zeeman 
splitting, and as in the perturbation approach, we keep 
only two states per rung : the singlet and the polarized 
triplet Then, we have computed up to range-8 

effective CORE interactions by solving exactly up to 8 
rungs. Despite having an infinite cluster expansion in 
Eq. {6]), previous studies have shown that, in many cases, 
the long-range effective interactions decay quickly so that 
they can be neglected beyond a certain range r. This is a 
necessary condition for any practical implementation of 
CORE and should be checked systematically. 

On Fig. [6l we plot the largest matrix element (in abso- 
lute value) of the range- r connected contribution h c r onn . 
For the coupling J x — 0.4, Jd = 0.2, we indeed observe 
a strong decrease of the amplitudes of the different pro- 
cesses as a function of the range of interaction r. This 
gives us confidence in the truncation beyond a certain 
range. However, the case of J x — 0.8 and Jd = 0.7 is 
an example where this CORE approach does not work, 
as will be discussed below. Here the matrix elements of 
jyconn remam substantial even for large ranges r. 

For the cases where the matrix elements of h c ° nn de- 
crease fast with increasing range r, we expect that CORE 
reproduces the low-energy physics of the system very 
well. In order to illustrate and strengthen this point, we 
have exactly solved different effective models obtained at 



a given truncation approximation. On Fig. [3 we show 
for the example of J x = 0.4, Jd — 0.2 that with increas- 
ing range, CORE results quickly converge to the exact 
ones. 

Since we have used the same basis for CORE as for 
perturbation theory in Sec. Ill Bt the range-2 effective 
Hamiltonian will also be a t — V model as in Eq. ([2]). 
But for CORE, the dependence of t, V, and /i on the 
different interactions is different and reads: 

t=^L V = Ess + p ± + ^±^ (7 ) 

M = -E S s - h - — 

Ess is the ground state energy of the 2x2 plaquette 
which is found to be: 



J± J x + Jd 
E SS = - — — 



(8) 



J\ +JZ + J}- J x Jd - J±(Jx + Jd) 



Consequently, the condition V/\t\ > 2 for the existence of 
a magnetization plateau at one half the saturation value 
translates into 



1 7/7 3 

< Jd/Jx < 



3 J x 



i + j x 



(9) 



This criterion is shown on Fig. \5\ together with the per- 
turbative result. For small interaction J x , one observes 
that the plateau phase boundaries given by Eqs. ^ and 
([9]) practically coincide; however, for increasing interac- 
tion J Xl i.e. when J x can no more be treated as a small 
perturbation, the two curves deviate from each other 
and CORE becomes more reliable. Note that a range-2 
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CORE calculation is quite simple, can be done analyti- 
cally by solving a 2 x 2 plaquette and already improves 
the accuracy with respect to perturbation theory. 

The limitations of this approach become evident for 
Jx ~ Jd > 0.7 where naive implementation of CORE 
fails to reproduce the absence of magnetization plateaux 
and the matrix elements of h c ° nn do not decrease with 
increasing range r, as shown for J x — 0.8, Jd — 0.7 in 
Fig. [6] There is a simple symmetry argument to under- 
stand why. For the specific case J x — Jd, the Hamil- 
tonian has many additional symmetries, namely the ex- 
change of the two spins on any given rung. Therefore, 
it can be shown that the effective hoppings are strictly 
zero at all orders in the CORE approach, which of course 
leads to a gapped insulating phase at half-filling. Indeed, 
in the CORE algorithm, it is necessary that the ground- 
state has a finite overlap in the reduced Hilbert space. 

In the vicinity of the line J x = Jd, this argument is 
no longer strictly valid but, for practical calculations, we 
observe that in the CORE calculation, many overlaps be- 
come very small, which results in effective models that 
are not accurate. Even if in principle, CORE could be- 
come accurate if longer-range interactions are taken into 
account, we think that in such cases, CORE loses its 
practical utility. 

As we will discuss in the following, when the chosen 
basis does not represent correctly the ground-state, prac- 
tical implementations of the CORE algorithm fail. When 
this is the case, several strategies are possible : (a) Keep 
other states for each block; (b) Keep more states for each 
block or (c) Change the block. We now turn to the dis- 
cussion of each case by focusing on what are the best 
block decomposition and/or states to keep. 



III. CHOICE OF BASIS 

CORE will be useful when one can keep a small num- 
ber of states per block and restrain to finite-range ef- 
fective interactions. However, checking the convergence 
is not always easy (except in ID as we have shown on 
Fig. [7]) so that it would be useful to have alternative in- 
formation. In Ref. [2(|, the authors proposed to used the 
reduced density matrix of the block as a tool to compute 
the relative weights of each block state. This procedure 
is similar to the density-matrix renormalization group 
method (DMRG)2^ and can both help to choose the cor- 
rect number of kept states, or indicate that a given block 
decomposition might not be appropriate. Moreover, this 
analysis can be done independtly of CORE since it relies 
on an exact calculation and is rather easy since it can be 
done on small clusters^. 



A. Rung density matrix 

A first illustration of this reduced density matrix 
weights is given on Fig. [5] We have considered a rung 
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FIG. 8: (Color online) Reduced density matrix weights for 
the vertical dimer block obtained on a 2 x 12 ladder as a 
function of J x and Jd along a path in the phase diagram. 
Calculations have been carried out for a magnetization of 1/2 
of the saturation value. 



embedded in a larger cluster and we trace out the other 
spins in the exact density matrix obtained with the 
m z = \m sat ground-state^ (GS) Y$>). By diagonalizing 
this reduced density matrix, we obtain the probablity of 
finding a certain block state, given that the overall system 
has a wavefunction For our choice of block, we ob- 
tain the weights of the rung states, namely the singlet \s) 
and the 3 triplets. We immediately observe that the sin- 
glet weight vanishes or is very small when J x ~ Jd > 0.7, 
which is precisely the region where a naive CORE cal- 
culation fails to reproduce the Haldane phase. In such 
cases, the exact GS has a vanishing or very small overlap 
in the naive CORE subspace, which explains the failure 
of our previous approach. 

On Fig. [5j a line indicates the region where the singlet 
weight becomes smaller than the S z = triplet weight. 
It corresponds precisely to the region with large J x ~ Jd 
where exact results have shown that there is no magne- 
tization plateau (see Fig. (2). In this region, one needs to 
use other strategies for CORE. 



B. Keeping 2 triplets 

In most of CORE approaches, the kept states have 
been taken as the lowest in energy on a single block. 
This choice is of course natural and is often the best one. 
However, in the case where both J x and Jd are close to 
J_l, we will argue that this is not the case. From the re- 
duced density matrix weight (Fig. [5]) , we have seen that 
the rung singlet state does not describe accurately the 
exact ground-state on a large system. Therefore, a first 
modification to the standard CORE algorithm would be 
to keep the size of the truncated basis to 2 but take the 
largest-weight states, namely |io) and \t+\) in this re- 
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FIG. 9: (Color online) Comparison of energy vs magnetiza- 
tion for J x — Jd = 1 on a 2 x 14 ladder given by ED and 
various CORE effective Hamiltonians obtained by choosing 
two triplets on each vertical rung and keeping up to range-r 
interactions. 



FIG. 10: (Color online) Comparison of energy vs magneti- 
zation for J x = 0.8 and Jd ~ 0.7 on a 2 x 14 ladder given 
by ED and various CORE effective Hamiltonians obtained by 
choosing one singlet and two triplets on each vertical rung 
and keeping up to range-r interactions. 



gion (the 3-fold degenerate triplets will have a Zeeman 
splitting in the presence of a magnetic field). As in the 
previous section, it is straightforward to compute ana- 
lytically the range-2 effective Hamiltonian that has again 
the form of t-V model for spinless fermions as in Eq. ([2|) 
with 



t 

V 



Jx + Jd 




J - +J 2 ri " J ' if J* ¥= Jd 



if J x — Jd 



(10) 



The crucial difference compared to the perturbative es- 
timate of Eq. is that the effective hopping of polar- 
ized triplets does not vanish anymore when J x — Jd- As 
a consequence, this effective mapping correctly predicts 
that V/t < 2 so that there is no magnetization plateau. 

Of course, we can go beyond range-2 approximation 
by including larger clusters in the CORE expansion. 
On Fig. [HI we consider the highly non-perturbative case 
where all couplings are equal and we observe a conver- 
gence of the CORE results towards ED data as we in- 
crease the range r. However, one needs to take into 
account longer-range effective interactions in the small- 
magnetization region and even for range-5 CORE the 
agreement with ED data is not satisfactory. 



difficulty lies respectively in (a) the very low density ma- 
trix weight of some states (b) keeping states that are not 
the lowest in energy. 

To resolve these difficulties, we decide to keep all the 
rung states except \t-i). Clearly, if we increase the size of 
the CORE subspace, we should get more accurate results 
but it will limit us in the possibility of studying such an 
effective model. On Fig. [TOl we have solved numerically 
various effective Hamiltonians obtained by keeping up to 
a given range r effective interactions. Although qualita- 
tively correct since we do not observe any finite plateau 
in this region, the CORE results converge slowly to the 
exact ones and give poor accuracy at small magnetiza- 
tion. 

Note that in the case J x — Jd > 0.8, the rung singlet 
weight vanishes (as shown in Fig. [H]) so that the CORE 
results are identical whether we keep two triplets and the 
singlet, or only two triplets. 

As a conclusion for the rung basis, for a given block 
decomposition, we propose that an efficient CORE imple- 
mentation should keep the first M low-energy states per 
block such that the total reduced density matrix weight 
is "large" . Of course, in some cases, we had to keep a 
rather large part of the total Hilbert space, which limits 
us both for numerical simulations and for analytic study 
of the effective model. Therefore, another route can be 
chosen by modifying the block decomposition. 



C. Keeping more states per block 

With the rung decomposition, we have observed that 
CORE is not very efficient if we restrict the truncated 
basis to singlet and polarized triplet. We have obtained 
better results by keeping the two dominant triplet states 
but the convergence was still poor. The origin of these 



D. Horizontal dimer block 

Since the rung basis does not give accurate results in 
the strongly frustrated regimes where all three couplings 
are of the same order, we decide to choose another block 
decomposition with horizontal dimers, keeping the sin- 
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FIG. 11: (Color online) Reduced density matrix weights for 
the horizontal dimer block obtained on a 2 x 12 ladder as 
a function of J x and Jd along a path in the phase diagram. 
Calculations have been carried out for a magnetization of 1/2 
of the saturation value. 



glet and polarized triplet: each block will again be rep- 
resented by a pseudo-spin 1/2. Fig. [TT1 illustrates the 
reduced density matrix weights of the horizontal dimer 
states^. 

By applying CORE, we obtain a new effective ladder 
with many-body effective interactions. In its cluster ex- 
pansion, the CORE Hamiltonian should in principle con- 
tain all kind of clusters interactions, including L-shape 
ones. However, because we have to deal with connected 
interactions, it can be shown that some cancellations oc- 
cur: for instance, since a given L-shape cluster appears in 
only one rectangular-shape cluster, its contribution ex- 
actly cancels out in the cluster expansion^!. We have 
carried out a CORE calculation including up to 6-rung 
interactions and we have solved these effective models 
exactly. 

On Fig. [T21 we compare range-4 and -6 calculations 
to exact results when J x = J<j = Jj_ = 1. Clearly, 
the horizontal dimer blocking scheme provides very 
good agreement with exact data. In particular, we 
have a good convergence of CORE data when including 
longer-range effective interactions. Moreover, although 
the rung basis wrongly predicted the existence of a 
magnetization plateau, here we observe a smooth energy 
vs magnetization curve and we have checked that there 
is no indication of any plateau, as is known from exact 
calculations in this region. 



At this stage, by choosing either a 2-site rung or leg 
blocking scheme, we are able to reproduce qualitatively 
the whole phase diagram for couplings J x and J<j varying 
from to Jj_. This gives us confidence that CORE can be 
used to reduce the complexity of any microscopic model 
and still gives a correct description of the properties in 
the presence of a magnetic field. For instance, the effec- 



FIG. 12: (Color online) Energy vs magnetization obtained on 
a 2 x 12 ladder with isotropic couplings J m = Jd = Jj_ = 1. 
CORE calculations are done with : vertical rungs blocks, 
keeping \s) and |£+i) and including all effective interactions 
(corresponds to infinite range); horizontal dimer blocks, keep- 
ing \s) and |t+i) and including range-4 and range-6 effective 
interactions. Exact results (ED) are shown for comparison. 



tive models that we have obtained can be solved exactly 
on clusters twice as large as for the original model. 

Because of the large flexibility in the choice of the 
block, we now turn to another decomposition of the lad- 
der system. 



E. Plaquette basis 

Another possible block decomposition of a ladder con- 
sists in 4-site plaquettes. We start the CORE algorithm 
by classifying its 16 states: two singlets, three triplets 
and one quintet. In the presence of a sufficiently strong 
magnetic field, only the polarized components will be rel- 
evant so that we restrict ourselves to 6 states : the fully 
polarized quintet \Q), the three fully polarized triplets 
| 2a), I 2b), and |2b), and the two singlet states \Sa), 
and | Sb)- In order to make connection with the previ- 
ous sections, we rewrite those states in terms of the rung 
basis ones: 

the polarized quintet state \S = 2, S z = 2): 

\Q) = \t +1 t+i) = \t+i) ® |t+i); E Q = -(J X + Jx + Jd) 

the polarized triplet states \S = 1,S Z = 1): 



\T A ) =-^(|* t+i)-|t+ito» 
\T B ) = ^=(\st +1 )-\t +1 s)) 



\Tr 



V2 



(\st +1 ) + \t +1 s)) 



E Ta 



(J± — Jx — Jd) 
(J± -J x + J d ) 
( J± +Jx~ Jd) 
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FIG. 13: (Color online) Energy of the plaquette states defined 
in Eq. as a function of the interaction along the ladder 
J x and the diagonal interaction Jd (for h = and Jj_ = 1). 

the singlet states 15 = 0, S z = 0): 

\S A ) = \ss); E SA = ~{J^ + J x + J d )-1 
\S B ) = -^(|i + it_i) + |*_ii+i>-|*o<o» 

E Sb = -\{JjL + J* + Jd)+1 (11) 

where 7 = J ' J\ + J| + Jj - Jj_ J x - Jj_ J d - J x J d 

Note that, concerning the energies, the contribution of 
the magnetic field h has been omitted. Fig. [TBI shows the 
behavior of the energies as a function of J x and Jd along 
a path in the phase diagram. 

We have previoulsy emphasized the use of the reduced 
density matrix weights in order to correctly choose which 
and how many block states should be kept for a given 
CORE calculation. On Fig. [14] (a), we plot the weights 
of our chosen plaquette states as a function of the ladder 
couplings (for each state, we compute its weight for all 
Sl ot and only plot the largest value); therefore, a small 
value indicates that a given plaquette state is not relevant 
to describe the exact GS for any 5*'°', i.e. for any mag- 
netic field. Using this information, we can even reduce 
further the CORE basis for some parameters : if some 
weights are tiny or even zero (typically smaller than 5%), 
we can reduce the CORE basis from 6 states to only 3 
(one in each ^-sector), and still have a good accuracy. 
Of course, reducing the CORE basis allows us to solve 
exactly the effective model on much larger system sizes 
(up to 2 x 36 for the original model) . We can also use the 
reduced density matrix weights to check whether the 6- 
state basis correctly reproduces exact GS properties: on 
Fig. Q~J] (b) is plotted the total weight of these six states 
(we have computed the sum of these six weights for all 
Sl ot and we only show the smallest value). Since this 



FIG. 14: (Color online) Reduced density matrix weights of 
the six chosen plaquette states (see Eq. 1 1 1 1> as a function of 
J x and Jd, obtained by ED from the exact ground-state of 
a 2 x 12 ladder with various S* *: (a) largest weight of the 
six plaquette states for all S' ot ; (b) minimum over S*°* of the 
cumulated weight of the six states. See text for details. 

total weight is larger than 60% in all the phase diagram 
(couplings and magnetic field), we expect that CORE 
effective interactions should decay quickly with distance. 

Then, by solving exactly the effective models, we can 
compare the energy vs S z curve to the exact one. On 
Fig. 1151 we plot our data obtained from various couplings 
throughout the phase diagram; CORE calculations are 
done either with 3 states (one per S z sector with the 
largest reduced density matrix weight) or all six states 
according to Fig.Q3] In all cases, the CORE convergence 
is very good and range-4 calculations are very accurate^. 
Being confident in our effective models, we can compute 
the magnetization curve on much larger systems and we 
show on Fig. [TB] typical plots showing either the presence 
or absence of a magnetization plateau at half-saturation, 
in full agreement with exact results. As a side remark, 
since our effective models are not necessarily particle- 
hole symmetric (in the bosonic language) , we observe on 
Fig. [TBT a) that the magnetization curve behaves differ- 
ently close to m z ~ and m z ~ m sat ; in particular, 
the precise shape close to the saturation value converges 
quite slowly with the range r of the effective CORE in- 
teractions. 

As a conclusion for the frustrated ladder, we have been 
able to obtain reliable effective Hamiltonians for various 
block decompositions (vertical or horizontal dimers, or 
plaquette). For some parameters, some choices are bet- 
ter than the others: for instance, in the non-perturbative 
frustrated regime where all three couplings are of the 
same magnitude, we have found that, with a plaquette 
decomposition and keeping only three states per block 
(not always the low-energy ones), we had a quick con- 
vergence of the CORE effective interactions that lead to 
a high accuracy, while CORE calculations are more dif- 
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FIG. 15: (Color online) Comparison of energy vs magneti- 
zation given by ED and CORE Hamiltonians keeping up to 
range-2 or 4 effective interactions for a 2 x 14 ladder with 
various couplings, (a-c) for CORE, we keep only 3 states 
per plaquette, which have the largest reduced density matrix 
weights (see Fig. I14J1 ; (d) for CORE, we keep six states per 
plaquette (see Eq. (fTTjl ). 
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FIG. 16: (Color online) Comparison of magnetization curve 
for 2xL ladder obtained with ED and range-3 CORE cal- 
culations with (a) three kept plaquette states for J x = 1 and 
Jd — where there is no plateau; (b) six kept plaquette states 
for J x — 0.7 and Jd = 0.5 where there is a finite plateau. 



ficult starting from rungs. The reduced density matrix 
criterion gives us two useful informations : (i) it gives a 
systematic way to locate the energy cut-off and fix how 
many low-lying states per block should be kept in the 
CORE approach; (ii) it can also indicate that relevant 
block states are not necessarily the low-energy ones. 



IV. TWO DIMENSIONAL SYSTEM 

A. Anisotropic case 

We now consider the 2-dimensional (2D) bilayer anti- 
ferromagnetic spin dimer XXZ model which is given by: 

Uxxz = ~hJ2 S *,i + ^X)S M .S 2li + (12) 

a, i i 
+ J (^o,t^Qj + ^a,i^a,j + ^Sa,iSa,j) J 

where a = 1,2 denotes the layer index. The interlayer 
coupling J is the largest one and has SU(2) symme- 
try, while the intra-layer coupling J is taken with an 
anisotropy A ( A = 1 corresponds to the isotropic SU (2) 
case). Such a model has been recently introduced and 
studied numerically with a QMC techniqu e) 39 ' 40 . As for 
the ladder model, it can be convenient to use the particle 
language, where the effective triplets behave as hardcore 
bosons that can have a solid (plateau region), superfluid 
(no plateau), or even a supersolid phase with both super- 
fluid and solid order parameters. QMC simulations have 
shown that this model exhibits all these phases^ 9 -, includ- 
ing a large half-saturated magnetization plateau region. 

Analogously to the perturbation theory that has been 
applied on the frustrated 2-leg ladder in section Iff Bl we 
can carry out a perturbation calculation, and then derive 
an effective 2D (hardcore) bosonic t — V model : 



H, 



tv 



(ij) (ij) » 



(13) 



or equivalently an effective XXZ spin-1/2 model. Again, 
we restrict ourselves to the singlet |s) and the polarized 
triplet |t+i) on each rung in order to describe the system 
in a magnetic field close to m = 1/2. 

One finds the following set of parameters: 

/' A /' 

t=- V=— n=J-h (14) 

This means that in perturbation theory, one finds a finite 
plateau (superfluid- insulator transition) when V/\t\ = 
A > A c = 2, which is independent of J and J' (for 
J' 7^ 0) . This line is drawn in Fig. [T7] and deviates from 
the QMC point. 

Performing the simplest CORE approach we use the 
same two states (|s) and |t+i)). Restricting to range-2 
interactions, the parameters for the effective t—V model 
can be easily derived: 



t 



.r_ 

2 



V = 

J 

~2 



AJ' 

— + 

- E G s 



E, 



GS 



3J 

T 



(15) 



Here Eqs is the ground state energy of the two dimers. 
We have also calculated the reduced density matrix 
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FIG. 17: (Color online) Critical anisotropy A c above which 
there is a plateau phase as a function of J' for perturbation 
theory and the simplest CORE approach for J = 1. The 
QMC point: J'/ J = 0.29 and A c = 3.2 is taken from Ref.H. 



weight for this case^ and observed that the singlet and 
the polarized triplet represent more than 90% of the to- 
tal weight in the region of < J' < 0.5 and < A < 4, 
which gives us confidence in the reliability of the effec- 
tive model. We also noted that the weights depend only 
marginally on J' and A (data not shown). For instance, 
the weights for A = 1 can be retraced in Fig. [TH] in the 
sector where J a = 0. 

Comparing with the perturbative result of Eq. (fT4)l . 
one sees that t remains unchanged whereas V and fi are 
modified. Consequently, the criterion for a superfluid- 
insulator transition (V/|t| = 2) is no longer independent 
of J and J': 



£ , 2E GS + 3J 

\t\- 2 - Ac+ J' 



(16) 



This expression only coincides with the result of the 
perturbative approach for J' — > 0. As is shown in 
Fig. |T7l the critical value of the anisotropy A c increases 
monotonously with J' if we set J = 1 . Above this curve, 
the CORE approach predicts a solid phase which means 
that the system exhibits a magnetization plateau at 1/2 
of its saturation value. For comparison, we have added 
the QMC result^ indicating a solid phase (i.e. finite 
plateau) when J' J J — 0.29 and A > 3.2, which is re- 
markably close to our simple estimate. Thus this CORE 
approach with two dimers which is modest and very sim- 
ple to carry out is a reliable tool to predict the occurence 
of a solid phase and gives a much better agreement with 
exact results than perturbation theory. 

Moreover, in the numerical study performed by Ng 
and Lee^, supersolid behaviour has been shown to oc- 
cur close to the solid phase. We believe that it would be 
desirable to derive effective bosonic models showing such 
a rich phase diagram, particularly for frustrated models 



which are not accessible by QMC simulations due to the 
negative-sign problem. Our simple t — V model of Eq. @ 
does not have a supersolid phase but exhibits phase sep- 
aration instead 41 . For that reason, we extend the range 
of our CORE approach to range-4, i.e. we consider a 
system of 2 x 2 dimers 4 ^. We then obtain an effective 
hardcore bosonic model containing all possible interac- 
tions allowed by symmetry (i.e. conserving the particle 
number) : 

C n 



+ \{b\bj+ +h.c.)+t?>(blb k + bjbj + h.c.) 
+ t^ibjblbibj + h.c.) + tf\b\b%b k + O +h.c.) 



"2 \ u l u k 1 J 
*( 3 )^t 



+ t\ >(b\bi(nj + n k - 2 nj n k )+ O +h.c.) 
+ tf\blb k ( nj +ni- 2n jni )+ (3) 
V (1) 

+ -^(n l n 1 +0) + V 2 2 (n l n k + n J n l ) 

r(l). 



+ V3- '(riinj(n k + ni - 2n k n t )+ O) 
t^\b\bjn k n 
V4mnjn k ni 



+ t[ 5) {b\bjn k ni+ O +h.c.)+t i ?\b\b k n j ni+ O) 



(17) 



Here we use the following notation: the sum X) ! n fc goes 

i j 

over the four plaquette sites. The circle O indicates 
that all contributions of one kind are included: E.g. 

Ejn* (4?(^+ O)) = l Mb\b 3 + b]b k + 6+6, + bjbi). 

In the region of interest, i.e. J' /J = 0.3 and < A < 

4, although all parameters are non-vanishing, we find 

that the dominant terms are the following: the nearest 

(1) (' 

and next-nearest neighbor hoppings t\ and t± 



(2) _ .(6) 



the nearest and next-nearest neighbor repulsion V 2 

(2) 

and Vn ■ The dependence of these parameters on the 
anisotropy can be seen in Fig. [TS] together with the 
range-2 parameters t and V that have been presented 
in Eq. (fT5")) . For the range-2 data we note that the tran- 
sition from the superfluid phase to the solid phase occurs 
at A ~ 2.8 for a inter-dimer coupling of J' = 0.3. 

Being able to compute a more refined effective Hamil- 
tonian, one could wonder about the possible effects of 
higher-order terms. This model is very complicated and 
there is no simple criterion to predict the existence or 
not of a plateau phase. Nevertheless, we believe that our 
range-4 effective Hamiltonian will also possess a plateau 
phase for large A. First of all, since this phase is gapped, 
we expect it to be robust to small additional interac- 
tions. Moreover, starting from a half-filled solid phase, 
the dominant longer-range terms have the following ef- 
fect: as shown in Fig. [TH1 the diagonal density interac- 
(2) 

tion V 2 < enhances the stability of the solid phase, 

while the diagonal hoppings and tf\ that would fa- 
vor a superfluid phase over the solid, remain small. This 
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FIG. 18: (Color online) For fixed J = 1 and J' = 0.3: t 
and V as function of anisotropy A for range-2 CORE (see 
Eq. (|15p ): dominant parameters: t' 1 ', V^ 1 ', and V^ 1 

also as function of A for range-4 CORE (see Eq. IfTT}). 



FIG. 19: (Color online) Reduced density matrix weights for 
the frustrated 2D Heisenberg bilayer obtained on a 2 x 16 
system as a function of J x and Jd along a path in the phase 
diagram. Calculations have been carried out for a magneti- 
zation of 1/2 of the saturation value. 



qualitative argument could be checked quantitatively by 
QMC simulations, for instance with a simpler effective 
model with no minus-sign problem. 

Moreover, the presence of diagonal single-particle hop- 

(2) 

ping t\ and other terms should also stabilize a super- 
solid phased. In particular, removing a particle from the 
half-filled solid creates a defect that can propagate due 
to diagonal hopping, which results in a superfluid order 
parameter coexisting with the solid, i.e. a supersolid. 
On the contrary, if we add a particle to the half-filled 
solid, our range-4 CORE effective Hamiltonian does not 
allow for diagonal hopping since the two processes de- 
scribed by 4 and 4 m Eq. (fT7|) exactly cancel out 
(tj 2 ' + = 0). From this observation, we expect that 
a supersolid phase is stable (respectively unstable) for 
filling smaller (respectively larger) than 1/2, which is in 
perfect agreement with the findings of Ref. [39|. 

It is also interesting to note that among the longer- 
range effective interactions, correlated-hopping terms 
could give rise to new phases^ and could be relevant for 
some geometries (e.g. Shastry-Sutherland lattice); here, 
we have found that their amplitudes remain very small. 



B. Frustrated case 

Let us now regard the 2D antiferromagnetic bilayer 
model without anisotropy A but with a frustration Jd 
between the two layers. This system corresponds to the 
2D generalization of the frustrated Heisenberg ladder of 
Eq. |T]). Fig. [T^l illustrates the reduced density matrix 
weights^ 

Carrying out a CORE calculation we can derive a,t — V 
model on two dimers (see Eq. (|13[1 ) or a more compli- 
cated hardcore boson model including up to four dimers 



interactions of the form of Eq. (fTT)) . 

In the range-4 CORE calculation, unlike in the pre- 
vious anisotropic case, we find that only the nearest- 
neighbor hopping 4 and repulsion V 2 f are relevant 
since the other terms are at least one order of magnitude 
smaller. These parameters are plotted in Fig. [201 Fur- 
thermore, one observes that for J x = 0.3 and < Jd < 
0.5 the range-2 calculation has already converged, i.e. 
the effective parameters are almost identical for range- 
2 and range-4 approximations. Considering the large 
reduced density matrix weights of the singlet and the 
polarized triplet for these parameters (the sum of their 
contributions exceeds 95% of the total weight, as shown 
in Fig. [19)1 . we expect that CORE converges fast in this 
region. 

Since the 2D or ID supcrfluid-insulator transitions oc- 
cur when V/\t\ — 2 which corresponds to a frustration 
of Jd ~ 0.1 for J x — 0.3, we expect a plateau phase re- 
gion similar to the one observed in the ladder model (see 
Fig. [5]). 

Depending on J x and Jd values, we find that higher- 
order terms may become important and could help sta- 
bilize a supersolid regime close to the plateau phase. A 
precise understanding of the effects of these many-body 
effective interaction is, however, beyond the scope of our 
study. 



V. CONCLUSION 

For frustrated ladders, naive perturbation in the rung 
basis completely fails to describe the disappearance of 
plateaux observed for large J x ~ Jd ~ J±- We have 
shown that, even with a simple range-2 CORE calcula- 



12 



1 


■ — 


~ - - t (range-2) 
- - V (range-2) 

■■■■ tj ' (range-4) 

-■■ V, (range-4) 
i 





0.2 0.4 0.6 0.8 



FIG. 20: (Color online) For fixed J x = 0.3 in the frustrated 
case: effective t and V parameters as a function of Jd for 
range-2 CORE (see Eq. ([3J ); and V 2 (1) also as function of 
J d for range-4 CORE (see Eq. (JT7J). 

tion keeping 2 states per rung, which can be done analyt- 
ically, CORE gives much more accurate boundaries for 
the plateau region, as seen when comparing exact data 
from Fig. [2] and CORE predictions on Fig. [5] For the 
CORE calculations, we have found that it can be crucial 
to use the information given by the exact reduced density 
matrix weights of the kept states in order to choose the 
best CORE basis, namely to answer the question: what 
is the best blocking scheme and how many block states 
should be kept ? Another advantage of the CORE calcu- 
lation is that its accuracy can be systematically improved 
by including longer-range effective interactions. By per- 
forming various block decompositions and keeping many- 
body effective interactions, we have shown that CORE is 
able to reproduce quantitatively the properties of frus- 



trated ladders in the presence of a magnetic field. The 
reduction of the Hilbert space allows us to solve exactly 
the effective model on much larger system sizes (up to 
2 x 36) compared to standard ED. 

We have also considered two-dimensional anisotropic 
or frustrated Heisenberg bilayers. In the non-frustrated 
anisotropic case, our CORE calculation improves over 
perturbation theory in locating the condition for a 
plateau formation and is compatible with a recent QMC 
study. Fine details such as the occurence of supersolidity 
can also be captured by computing longer-range effective 
interactions, such as diagonal hopping. In particular, we 
find that assisted hopping is not a relevant interaction 
for that geometry and these parameters. The frustrated 
bilayer was shown to have a similar effective model with 
similar amplitudes. Therefore, we predict that it will 
have a similar phase diagram as the anisotropic bilayer, 
containing superfluid, solid and supersolid regions. The 
advantage of such a model is that it respects SU (2) sym- 
metry and could be more realistic for materials descrip- 
tion. Physically, in the anisotropic case, the effective 
repulsion between hardcore bosons increases with grow- 
ing anisotropy, leading to an insulating phase, whereas in 
the frustrated case the effective hopping is also strongly 
reduced with growing frustration. 
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